UK Patent Application <«,GB o„2356455 „3,A 



(43) Date of A Publication 23.05.2001 



(21) 


Application No 9926929.2 


(51) 


- 

INT CL 






GOl V 1/28 


(22) 


Date of Filing 16.11.1999 








(52) 


UK CL (Edition S ) 






GIG GEL 


(71) 


Applicant(s) 








ueco-KraKia (Uixi umitaa 


(56) 


Documents Cited 




(tncorporateo in the United Kingdom) 




US 5657294 A US 5142501 A US 5060204 A 




Schlumberger House, Buckinghain Gate, GATWICK, 




US 5060203 A 




West Sussex, RH6 ONZ, Untted Kingdom 










(58) 


Field of Search 


(72) 


Inventor(s) 




UK CL (Edition R ) G1GGEL 




Richard Bale 




INTCL^ GOl V 1/28 




Gabrieia Dumitru 




ONUNE: EPODOC, WPI, JAPIO 


(74) 


Agent and/or Address for Service 








Marks & Clerk 








4220 Nash Court Oxford Business Park South, 








OXFORD, OX4 2RU, United Kingdom 







(54) Abstract Trtle 

Determination of fast and slow shear wave polarisation directions 

(57) A method of determining the polarisation directions of the fast and slow shear waves arising from 
shear wave splitting due to anisotropy, said directions defining a natural coordinate system, the method 
comprising the steps ot 

a) recording at least two components of each shear wave, in a recording coordinate system, 

b) calculating the value of e, being the angle of rotation between the natural coordinate system and the 
recording coordinate system, for which the Lp norm of the rotated traces is minimised if p is less than 2, 
or maximised if p is greater than 2. 
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Determinatioii of the Fast and Slow Shear Wave Polarisation Directions 

The invention relates to the detennination of the polarisation directions for the fast and 
slow shear waves arising from shear wave splitting due to anisotropy . 

5 

There are generally two types of seismic waves used in seismology, namely so-called 
"P-waves" or compressional waves in which the vibrations occur in the direction of 
propagation of the waves, and so-called "S-waves" or shear waves in which the 
vibrations occur in a direction generally orthogonal to the direction of propagation of 
10 the waves. 

A multicomponent geophone is a directional detector for seismic waves, which includes 
a vector measurement of the incoming wave. In the applications considered here, two of 
the geophone components are assumed to be aligned along arbitrarily chosen X and X 
IS directions, generally parallel to the surface of the earth. It is also assumed that the 
incoming shear waves will generally arrive vertically (i.e. perpendicular to the surface 
of the earth) from below the geophones. As a result the particle motion within the wave 
is generally parallel to the surface of the earth, and is detected by the X and Y geophone 
components. 

20 

Furtheraiore, as will be explained below, the incoming shear waves may contain two 
components which are polarised (in terms of the direction of vibration) in two 
orthogonal directions, SI (i.e. the fast shear Si propagation direction) and S2 (i.e. the 
slow shear S2 propagation direction), and which are separated from each other by a time 
25 delay. This specification is concemed with the determination of these two directions and 
the travel time delay between the corresponding shear waves. 

From ocean bottom or land multicomponent surveys using a P-wave source, it is 
possible to obtain measurements of the shear waves converted in the earth. These shear 
30 waves appear predominantly on the horizontal Q[ and Y) components of the 
multicomponent geophones. If the earth is isotropic with respect to the horizontal 
direction of wave motion, then a single shear arrival may be expected for each reflecting 
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interface. If however, as is often the case, the earth behaves anisotropically v^ith respect 
to the horizontal direction (for example, because a geological layer is polarised in a 
particular direction due to fracturing), then we can expect to record two separate shear 
wave arrivals from each reflecting interface, arriving at different times, having 
5 propagated with different velocities. These are usually termed the fast (Si) and the slow 
(S2) shear waves, corresponding to the first and the second amvals, respectively. They 
are also characterised by having differmt polarization directions (i.e. directions of 
particle motion in the horizontal plane), which in most cases are considered to be 
approximately orthogonal to each other. It is assumed that this is the case here. 

10 

The shear wave splitting phenomenon is illustrated in Figure 1, which depicts a shear 
wave arrival (S) that, at the start (A) of an anisotropic medium, splits into two separate 
shear waves (Si and S2)9 having different polarisation directions and propagating 
separately with differing velocities imtil the end (B) of the medium. If from (B) onwards 
15 the medium is supposed to be isotropic, the two polarised waves will continue to travel 
separately but with the same velocity until they impinge upon the recording geophones. 
The amplitudes recorded on each of the horizontal components of the geophone depend 
upon the orientations of the Si and 52 directions relative to the X and Y directions. 

20 Figure 1 gives a simple graphical description of the principle of shear wave 
birefringence, by only considering one anisotropic layer imbedded in an isotropic 
medium. However, in reahty there are many reflecting boundaries that give rise to a 
number of shear arrivals polarised in the SI and S2 directions. In addition, these SI and 
S2 directions can change between the different anisotropic layers. In the applications 

25 considered here, the SI and §2 polarisation directions are assumed to be constant with 
depth, over the analysing time window. 

Accordihig to the invention, there is provided a method of deteraiining the polarisation 
directions of the fast and slow shear waves arising from shear wave splitting due to 
anisotropy, said directions defining a natural coordinate system, the method comprising 
30 the steps of: 
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a) recording at least two components of each shear wave, in a recording 
coordinate system, 

b) calculating the value of 6 , being the angle of rotation between the natural 
coordinate system and the recording coordinate system, for which the Lp norm is 

5 minimised if p is less than 2, or maximised if p is greater than 2. 

In one embodiment of the invention, said value of 9 is determined by calculating the 
value of the Lp norai over a range of incrementally varying values of 9, and selecting 
that value of 9 for which the Lp norm is appropriately minimised or maximised. 

10 

In a further embodiment of the invention, p is 4, and the value of 9 is determined 
analytically from an equation derived by difTerentiating the Lp norm with respect to 9. 

The two recorded components of each shear wave are sampled, for example, at about 
IS 4ms intervals. 

Preferably, the fast and slow shear waves are recorded using two ordiogonal geophones, 
arranged generally parallel to the surface of the earth. 

The fast and slow shear waves may be produced fix)m a single source. 

20 Said source may be a P-wave source or it may be a single shear source. 

Said shear wave components are conveniently horizontal components. 

The invention also includes apparatus for carrying out the above method, and a 
25 computer readable mediimi carrjdng a computer program for carrying out the above data 
processing steps. 

The invention will now be more particularly described, by way of example only, with 
reference to the accompanying drawings, in which: 

30 
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Figure 1 shows shear wave splitting through an anisotropic medium; 

Figure 2a shows S\ and S2 traces (synthetically generated, as they would be recorded by 
5 geophones aligned with the natural coordinate system defined by the S1-S2 directions), 
with 30ms travel time delay between the fast and the slow shear waves, and in which §1 
is 30 degrees from the horizontal X axis of the recording system; 

Figure 2b shows X and Y recorded traces (synthetically generated, as they would be 
10 recorded by geophones aligned with the actual recording system, and corresponding to 
the measurements of arrivals in Figure 2a); 

Figure 3 shows pairs of and rotated traces after rotation of the X and Y recorded 
traces by angles ranging from 10 to 180 degrees, in which the Si S2 response is fully 
1 5 recovered when the angle used for the rotation is correct, in this case 30 degrees; 

Figure 4 illustrates that rotation of the X and Y axes to the Si and S2 directions is 
equivalent to moving trace samples, indicated by the stars, along circular paths with 
fixed distance from the origin; 

20 

Figure 5 shows a graphical comparison of constant LI, L2 and L4 norm contours; 

Figure 6 shows the result of applying the Li norm to the X and Y traces of Figure 2b for 
different rotation angles; it shows that the Li norm is minimised at 30 and 120 degrees, 
corresponding to the Si and S2 directions, respectively; 

25 

Figure 7 shows the result of applying the L4 norm to the X and Y traces of Figure 2b for 
different rotation angles; it shows that the L4 norm is maximised at 30 and 120 degrees, 
corresponding to the Si and S2 directions, respectively; 
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Figure 8 shows the estimation of the travel time delay by cross-correlation of with 
X^; the peak occurs at -SOms, meaning that the rotated component is in fhe fast SI 
direction, and that the lag between fhe Si and S2 shear waves is 3Qms. 

5 A number of methods have been proposed to determine the orientations of S| and S2 
shear waves and the corresponding travel time delay between them by using the 
recorded arrivals. Most of these methods rely upon having two shear sources with 
differing alignments, giving four independent measurements from two receiver 
components (e.g. Alford rotation, cf. Alford, 1986), The present specification is 
10 concerned with the case of a single source, typically generating P-waves, but possibly 
generating S waves. The SI and §2 directions are estimated using just two independent 
measurements from the horizontal X and X components of a multicomponent geophone. 

The problem is illustrated in Figures 2a and 2b* Figure 2a shows a simple example of Si 
15 and S2 arrivals in the **unknown'* natural SI -§2 coordinate frame, determined by the 
anisotropic layer from which reflection takes place. The Si and S2 traces are generated 
synthetically, as they would be recorded by geophones aligned with the SI and S2 
directions. Figure 2b shows the traces corresponding to the measurement of the arrivals 
in Figure 2a, recorded by geophones aligned with the actual X and Y recording 
20 directions. Following a convention in this field, the positive signals are "filled in" in 
black ink, and the negative signals are not In Figure 2a, each pair of pulses (one on S\ 
and the other on 82) corresponds to a reflection from a different anisotropic layer. The 
two pulses are separated by a time delay because, due to their different polarisation 
directions, they travel at different speeds within the anisotropic layer. 

25 

It can be seen that traces in the S1-S2 coordinate system is a '"simpler" more 
"parsimonious" representation of the particle motion than those in the X-X coordinate 
system, in the sense that each reflector only gives rise to a single event on the S| trace 
and a single event on the S2 trace, whereas there are two events on the X and Y recorded 
30 traces for each reflection (due to the fact that each one of the two shear waves are 
recorded by both the X and Y geophones). In seismic, and other signal processing 
fields, the simplest representation is sometimes referred to as the ''minimiun entropy" 
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rq>resentation, by analogy with thennodynamics. For this reason, the method described 
here may be termed the Minimum Entropy Rotation (MER). 

Parsimony can be measured by using a "norm". As an example, let us consider the Li 
5 norm. For each trace (typically sampled at about 4ms intervals), the Li norm is 

computed by summing the absolute values of die trace samples. The most parsimonious 
form of the trace is that which has the lowest Li norm. This principal is used in seismic 
analysis to perform "sparse spike inversion". 

10 The Li norm example can be generalised by using the linear p-norm Lp which, for a 
trace X ={xi, X2, Xn}, can be written as: 



where in all the following formulas 7V^ is the number of samples in the trace, and p is a 
real number. In the case of 2-component data, with traces X and Y, the norm can be 
1 5 calculated using samples from both traces: 



Another measure of parsimony previously used in seismic (for single component data 
only) is the *Varimax norm" (Wiggins, 1977) defined by the ratio of the fourth moment 
to the second moment squared: 




(1) 




(2) 



N 



20 V(X) = 




(3) 
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In statistics this is referred to as **Kurtosis". In terms of linear norms, the varimax norm 
is the fourth power of L4/L2. The most parsimonious result is obtained when the 
varimax norm attains its maximimi value, 

S Wiggins first used this measure of parsimony to determine the parameters of the 
deconvolution operator that best improves trace resolution. In this view, a seismic 
deconvolution operator is determined such that, when applied to a seismic trace, it 
produces an output with the greatest varimax norm. This method is known as the 
"minimum entropy deconvolution" (Wiggins, 1977). 

10 

In the following, we describe and compare different approaches developed to determine 
the shear wave splitting parameters for a computer generated model of the Si and S2 
traces. We start with a scanning approach over the recorded traces rotated by differing 
angles, and then make use of the Lp norm defined above, for different values of p. 

15 

Figure 3 shows the result of rotating the X and Y traces of Figure 2b by differing 
angles. It can be seen that the traces corresponding to the 30 degrees rotation are similar 
to those shown in Figure 2a, Thus, when the correct rotation angle is applied, the Si and 
S2 response is recovered. 

20 

For estunation of an imknown pair of axes SI and SIj which in most cases are rotated 
relative to the X-Y frame, we calculate the combined norm of the X and Y traces. We 
may make use of any norm other than L2, as this is invariant under rotation, and thus not 
useful. If using a Lp norm with p<2, it is necessary to find the angle 0 that minimises the 
25 norm, while if a Lp norm with p>2 is used, it is necessary to maximise the norm in order 
to find the angle 0 corresponding to the desired rotation, namely the angle 9 between the 
Si -S2 and X-Y coordinate systems. 

The principle is illustrated using the Li norm. 

30 The Li norm of the i-th trace sample of the rotated traces, X^ and Y^, is written for each 
angle 0 as follows: 
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/O) . (e ) = (e )| + ^yf (9 1 = cosO + sinG | + {x^ sinG - 3;,. cosG | , 



(4) 



10 



15 



20 



where x\ and y,- are the recorded trace samples of the X and Y traces respectively, G is 
the proposed angle between X-Y and Si-S2 coordinate systems, and xj^ and yj^ are the 
rotated trace samples. 

The total Li norm over a time window with samples is then written as: 



This norm is minimised when the rotated traces are in their simplest form, which is 
when the angle G is that of the SI or S2 directions, relative to the X or Y axes of the 
recording system* This can be seen by reference to Figure 3. 

Figure 4 depicts the dependence on the rotation angle of the Li norm |xj^|+|yi^|. The 
stars represent trace samples from the Si and S2 polarised arrivals* Applying a rotation 
of the X and Y axes to SI and S2 directions is equivalent to moving the trace samples 
around circles towards the SI and S2 axes. As the samples are rotated, their distance 
from the origin (x^f+{y\^f remains fixed. However, the sum of their absolute values 
varies with the rotating angle, attaining its minimum value after rotation by G, the angle 
between the X axis and the Si axis, or by 6+90, the angle between the X axis and the S2 
axis. This is used to estimate the directions of the SI and S2 axes corresponding to the 
fast Sj and slow S2 shear waves, respectively. 

The Li contours of constant |xi VlYi^l form diamond shapes as shown in Figure 4. The 
minimum Li value for a fixed distance from the origin occurs when the rotation is equal 
to G. Note that the L2 contour of constant (Xi^)^+(y,*^)^ is in fact the confining circle. 




(5) 



The method described above exploits the fact that rotation moves data samples along 
circles (i.e. constant distance from the origin). Therefore, for the L2 norm, the traces 
rotated to the differing angles are characterised by the same norm value, and so this 
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norm cannot be use to detemiine the shear wave splitting parameters. The Li nomi 
contours are diamond shaped, with their comers on the X and Y axes, attaining their 
minimum value (for a fixed distance from the origin) when data is rotated by the angle 
between X and SI or S2. This is generally true for all the Lp norms, having p<2. On the 
5 other hand, for p>2 the contours are more square shaped with the flatter sides on the X 
and Y axes. These nonns attain their maximum value when data is rotated by the angle 
between X and SI or §2. The Lp norm shapes are illustrated for p=l,2,4 in Figure 5. 

Figure 6 shows the value of the total Li norm, as given by equation (5) and applied to 
10 the X and Y traces of Figure 2b, plotted against the rotation angle. This displays two 
clear minima at the angle values of 30 and 120 degrees, corresponding to SI and S2 
directions. 

In practice, the L] norm is not die most convenient one to use, as finding the solution of 
15 the rotation angle requires the bmte force scanning approach described above. That is, it 
is necessary to calculate the norm for each increment of, say 1 degrees, in order to find 
the value of 6 corresponding to the mirumum value of the norm. A better choice is the 
L4 norm, for which an analytical treatment is possible. 

20 The L4 norm of the i-th trace samples of the rotated traces, X^ and Y^, is written for 
each angle 6 as follows: 

/ W,(e ) = y (e )f + (e )f = [x, cose + sine f + ^ sine - cose f (6) 
The total L4 norai over a time window with samples is then given by: 

^4(e)=i:/H(e) (7) 

1=1 

25 This can be differentiated with respect to the 0 angle to find that the L4 norm has its 
extreme values for 



9 =±arcsin 



(8) 



where: 




(9) 



1=1 



There are eight solutions given by equation (8), four of which are spurious. Of the 
5 remaining four, two give the L4 minima and two give the L4 maxima. In order to 
identify the desired solution we substitute them into the L4 norm equation (7) and select 
one of the two valid solutions that produces the same maximum value for L4. These two 
solutions represent the angle between the X axis and the Si and S2 directions. It is not 
important which of the two solutions are selected at this stage, as the next step (i.e. 
10 cross-correlation) will help to distinguish between the fast and slow shear directions. 

Figure 7 shows the value of the total L4 norm, as given by equation (7) and applied to 
the X and Y traces of Figure 2b, plotted against the rotation angle. This displays 
maxima at the angle values of 30 and 120 degrees, corresponding to SI and S2 
15 directions. 

There are some important advantages in using the Minimum Entropy Rotation method 
described l;ere. For single source component data, most of the existing methods used for 
estimating the shear wave splitting parameters assume that the shear wave is equally 

20 reflected and transmitted in both the SI and S2 directions. In general, this assumption is 
not always true, hence the strength of the invention presented here which does not 
require making this assumption. It is our observation that the method is also not very 
sensitive to noise. This could be explained as follows: due to additive noise, the samples 
.may not lie on the SI and S2 axes, but can be expected to have a mean value which 

25 approximately does; provided enough samples are included in the analysis and the noise 
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is not correlated with a particular polarisation direction, the estimation will still behave 
well 

As well as the angles giving the polarisation directions of S| and S2, an important shear 
5 wave splitting parameter is the travel time delay between the fast and slow shear waves. 
This may be obtained by cross-correlation of the two rotated traces at one of the 
minimum Li or maximum L4 positions. The maximum cross-correlation output is 
picked to determine: 

1. Which of the two traces contains the fast shear (i.e. the first arrival). 
10 2. The travel time delay between the fast and slow shear arrivals. 

Figure 8 shows the cross correlation of the with the traces rotated to 30 degrees. 
The negative time of the peak indicates that the trace is the fastest, whilst the time of 
30ms indicates the delay time between the fast and the slow shear waves. There is a 
1 5 secondary peak, which is due to a second event on the input traces. 
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CLAIMS 

A method of determining the polarisation directions of the fast and slow shear 
waves arising from shear wave splitting due to anisotropy, said directions 
defining a natural coordinate system, the method comprising the steps of: 

a) recording at least two components of each shear wave, in a recording 
coordinate system, 

b) calculating the value of 9 , being the angle of rotation between the natural 
coordinate system and the recording coordinate system, for which the Lp norai is 
minimised if p is less than 2, or maximised if p is greater than 2. 

A method as claimed in claim 1, wherein said value of 6 is determined by 
calculating the value of the Lp norni over a range of incrementally varying values 
of 9, and selecting that value of 9 for which the Lp nomi is appropriately 
minimised or maximised. 

A method as claimed in claim 1, wherein p is 4, and wherein the value of 9 is 
determined analytically from an equation derived by differentiating the Lp norm 
with respect to 9. 

A method as claimed in any preceding claim, wherein the two recorded 
components of each shear wave are sampled at about 4ms intervals. 

A method as claimed in any preceding claim, wherein the fast and slow shear 
waves are recorded using two orthogonal geophones arranged generally parallel 
to the surface of the earth. 

A method as claimed in any preceding claim, which involves producing the fast 
and slow shear waves from a single source. 
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A method as claimed in claim 6, wherein said source is a P-wave source. 

A method as claimed in any preceding claim, wherein said components are 
horizontal components. 

A method of determining the polarisation directions of the fast and the slow 
shear waves, said method being substantially as hereinbefore described with 
reference to the accompanying drawings. 

A shear wave analysis apparatus arranged to carry out the method of any 
preceding claim, the apparatus comprising shear wave detection means for 
recording said at least two components of each shear wave, and data processing 
means for calculating said value of 6. 



A computer readable medium carrying a computer program adapted to cause a 
computer to carry out the data processing steps of any of claims 1 to 9. 



